#include "RICHPhysicsList.hh"
//ConstructParticle()
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
//ConstructProcess()
#include "G4ProcessManager.hh"
// ConstructDecay()
#include "G4Decay.hh"
//ConstrucEM()
#include "G4Cerenkov.hh"
#include "G4Scintillation.hh"
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"
#include "G4eMultipleScattering.hh"
#include "G4hIonisation.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4MuMultipleScattering.hh"
#include "G4hMultipleScattering.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4OpAbsorption.hh"
#include "G4OpRayleigh.hh"
#include "G4OpBoundaryProcess.hh"
//Birks Correction to Scintillation
#include "G4LossTableManager.hh"
#include "G4EmSaturation.hh"

RICHPhysicsList::RICHPhysicsList() :
	QGSP_BERT() {
	;
}

RICHPhysicsList::~RICHPhysicsList() {
	;
}

void RICHPhysicsList::ConstructParticle() {
	QGSP_BERT::ConstructParticle();
	//	ConstructBosons();
	//	ConstructLeptons();
	//	ConstructMesons();
	//	ConstructBaryons();
}

void RICHPhysicsList::ConstructBosons() {
	G4Gamma::GammaDefinition();
	G4OpticalPhoton::OpticalPhotonDefinition();
}

void RICHPhysicsList::ConstructLeptons() {
	G4Electron::ElectronDefinition();
	G4Positron::PositronDefinition();
	G4NeutrinoE::NeutrinoEDefinition();
	G4AntiNeutrinoE::AntiNeutrinoEDefinition();
	G4MuonPlus::MuonPlusDefinition();
	G4MuonMinus::MuonMinusDefinition();
	G4NeutrinoMu::NeutrinoMuDefinition();
	G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
}

void RICHPhysicsList::ConstructMesons() {
	G4PionPlus::PionPlusDefinition();
	G4PionMinus::PionMinusDefinition();
	G4PionZero::PionZeroDefinition();
	G4KaonPlus::KaonPlusDefinition();
	G4KaonMinus::KaonMinusDefinition();
	G4KaonZero::KaonZeroDefinition();
	G4KaonZeroShort::KaonZeroShortDefinition();
	G4KaonZeroLong::KaonZeroLongDefinition();
}

void RICHPhysicsList::ConstructBaryons() {
	G4Proton::ProtonDefinition();
	G4AntiProton::AntiProtonDefinition();
	G4Neutron::NeutronDefinition();
	G4AntiNeutron::AntiNeutronDefinition();
}

void RICHPhysicsList::ConstructProcess() {
	QGSP_BERT::ConstructProcess();
	//	AddTransportation();
	//	ConstructEM();
	ConstructOp();
}

// Test if Particle Decays and Define Decay Process
void RICHPhysicsList::ConstructDecay() {
	G4Decay *theDecayProcess = new G4Decay();
	theParticleIterator->reset();
	while ((*theParticleIterator)()) {
		G4ParticleDefinition *particle = theParticleIterator->value();
		G4ProcessManager *pmanager = particle->GetProcessManager();
		if (theDecayProcess->IsApplicable(*particle)) {
			pmanager->AddProcess(theDecayProcess);
			pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
			pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
		}
	}
}

//Define Processes for non-Optical Particles
void RICHPhysicsList::ConstructEM() {
	theParticleIterator->reset();
	while ((*theParticleIterator)()) {
		G4ParticleDefinition *particle = theParticleIterator->value();
		G4ProcessManager *pmanager = particle->GetProcessManager();
		G4String particleName = particle->GetParticleName();

		//Construct Gamma Processes
		if (particleName == "gamma") {
			pmanager->AddDiscreteProcess(new G4GammaConversion());
			pmanager->AddDiscreteProcess(new G4ComptonScattering());
			pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());
		}
		//Construct Electron Processes
		else if (particleName == "e-") {
			pmanager->AddProcess(new G4eMultipleScattering(), -1, 1, 1);
			pmanager->AddProcess(new G4eIonisation(), -1, 2, 2);
			pmanager->AddProcess(new G4eBremsstrahlung(), -1, 3, 3);
		} else if (particleName == "e+") {
			pmanager->AddProcess(new G4eMultipleScattering(), -1, 1, 1);
			pmanager->AddProcess(new G4eIonisation(), -1, 2, 2);
			pmanager->AddProcess(new G4eBremsstrahlung(), -1, 3, 3);
			pmanager->AddProcess(new G4eplusAnnihilation(), 0, -1, 4);
		} else if (particleName == "mu+" || particleName == "mu-") {
			pmanager->AddProcess(new G4MuMultipleScattering(), -1, 1, 1);
			pmanager->AddProcess(new G4MuIonisation(), -1, 2, 2);
			pmanager->AddProcess(new G4MuPairProduction(), -1, 4, 4);
			pmanager->AddProcess(new G4MuBremsstrahlung(), -1, 3, 3);
		}
	}
}
//Define Processes for Optical Photons
void RICHPhysicsList::ConstructOp() {
	G4Cerenkov *cerenkovProcess = new G4Cerenkov("Cerenkov");
	G4Scintillation *scintillationProcess = new G4Scintillation("Scintillation");
	G4OpAbsorption *absorptionProcess = new G4OpAbsorption();
	G4OpRayleigh *rayleighProcess = new G4OpRayleigh();
	G4OpBoundaryProcess *boundProcess = new G4OpBoundaryProcess();

	//Cerenkov Processs
	//	cerenkovProcess->SetMaxNumPhotonsPerStep(1);
	//	cerenkovProcess->SetMaxBetaChangePerStep(.001);
	//	cerenkovProcess->SetTrackSecondariesFirst(false);

	//Scintillation
	scintillationProcess->SetScintillationYieldFactor(1.);
	scintillationProcess->SetTrackSecondariesFirst(true);

	//Birks Correction to Scintillation Process
	G4EmSaturation *emSaturation = G4LossTableManager::Instance()->EmSaturation();
	scintillationProcess->AddSaturation(emSaturation);

	//Boundary Process
//	G4OpticalSurfaceModel themodel = unified;
//	boundProcess->SetModel(themodel);

	theParticleIterator->reset();
	while ((*theParticleIterator)()) {
		G4ParticleDefinition *particle = theParticleIterator->value();
		G4ProcessManager *pmanager = particle->GetProcessManager();
		G4String particleName = particle->GetParticleName();
		if (cerenkovProcess->IsApplicable(*particle)) {
			pmanager->AddProcess(cerenkovProcess);
			pmanager->SetProcessOrdering(cerenkovProcess, idxPostStep);
		}
		if (scintillationProcess->IsApplicable(*particle)) {
			pmanager->AddProcess(scintillationProcess);
			pmanager->SetProcessOrderingToLast(scintillationProcess, idxAtRest);
			pmanager->SetProcessOrderingToLast(scintillationProcess, idxPostStep);
		}
		if (particleName == "opticalphoton") {
			pmanager->AddDiscreteProcess(absorptionProcess);
			pmanager->AddDiscreteProcess(rayleighProcess);
			pmanager->AddDiscreteProcess(boundProcess);
		}
	}
}

void RICHPhysicsList::SetCuts() {
	SetCutsWithDefault();
}
